
%%%READING FROM PARAMETERS.TXT
fpar=fopen('parameters.txt', 'r');
Parm=textscan(fpar, '%c %s %c %f%*[^\n]',22, 'headerLines', 1);
Parm{4};
fclose(fpar);
fpar2=fopen('parameters.txt', 'r');
Data=textscan(fpar2, '%s %f%*[^\n]',24, 'headerLines', 24);
Data{2};
fclose(fpar2);

ndt0=Data{2}(2)
ndg0=Parm{4}(14)
np0=Parm{4}(15)
sample_ndg=Parm{4}(18)
sample_ndt=Parm{4}(19)
deltat=Data{2}(7)
deltag=Data{2}(6)
sample_inc_t=Data{2}(16)
sample_inc_g=Data{2}(17)
simulation_steps_t=Data{2}(22)
simulation_steps_g=Data{2}(23)
sampled_particles=Parm{4}(20)
simulation_steps_p=Data{2}(24) %real simulated particles. must be equal to sampled_particles
w0=Parm{4}(22)

lightC_au=137.037;

fid=fopen('output1.txt', 'r');
subplot(2,3,1);
sp2=textscan(fid, '%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f','delimiter', ' ', 'headerLines', 6);
plot(sp2{1}, sp2{4},'LineWidth',1)

%axis([0 4*1024 -0.5 0.5])
xlabel('Time')
title('Transmitted field Ey')
Y2 = fft(sp2{4});
length(Y2)
%dt=sample_inc_t;
%ndt=simulation_steps_t;


dw=2*pi/ndt0/deltat
wmax=pi/deltat
subplot(2,3,2);
w=0:dw:wmax-4*dw;
length(w/w0)
length(log10((abs(Y2(1:length(Y2)/2)))  ))
plot(w/w0,log10((abs(Y2(1:length(Y2)/2)))  ))

%Pyy = Y.* conj(Y);
%plot(s{1},Pyy)
xlim([0 20])
ylim([0 max(log10((abs(Y2(1:length(Y2)/2)))))])
grid
title('Spectrum of the tansmitted field')
xlabel('frequency in harmonics')

fclose(fid);



subplot(2,3,4);
plot(sp2{1}, sp2{3},'r','LineWidth',1)

%axis([0 4*1024 -0.5 0.5])
xlabel('Time')
title('Reflected field Ey')

Y1 = fft(sp2{3});

%w0=5.700000e-02;

%dw=2*pi/ndt/dt;
%wmax=pi/dt;
subplot(2,3,5);
%w=0:dw:wmax-3*dw;
plot(w/w0,log10((abs(Y1(1:length(Y1)/2)))), 'r')
xlim([0 20])
ylim([0 max(log10((abs(Y1(1:length(Y1)/2)))))])
grid
title('Spectrum of the reflected field')
xlabel('frequency in harmonics')



subplot(2,3,3);
plot(w/w0,log10((abs(Y2(1:length(Y2)/2))).^2))
title('Spectrum of the tansmitted field (in E^2)')
xlim([0 20])
ylim([0 max(log10((abs(Y2(1:length(Y2)/2))).^2))])
grid
xlabel('frequency in harmonics')

subplot(2,3,6);
plot(w/w0,log10((abs(Y1(1:length(Y1)/2))).^2), 'r')
title('Spectrum of the reflected field (in E^2)')
xlim([0 20])
ylim([0 max(log10((abs(Y1(1:length(Y1)/2))).^2))])
grid
xlabel('frequency in harmonics')
% figure
% subplot(2,2,1);
% plot(sp2{1}, sp2{6},'r','LineWidth',1)
% 
% xlabel('Time')
% title('Exleft')
% 
% subplot(2,2,3);
% plot(sp2{1}, sp2{7},'r','LineWidth',1)
% 
% xlabel('Time')
% title('Exright')
% 
% subplot(2,2,2);
% plot(sp2{1}, sp2{8},'g','LineWidth',1)
% 
% xlabel('Time')
% title('Ro0left')
% 
% subplot(2,2,4);
% plot(sp2{1}, sp2{8},'g','LineWidth',1)
% 
% xlabel('Time')
% title('Ro0right')

figure
subplot(1,2,1);
plot(sp2{1}, sp2{5},'r','LineWidth',1)
xlabel('Time')
title('Incident field')
subplot(1,2,2);
%px=0:1:10;
%py=zeros(1,10);
axis off
text(0.5,0.5,[' E_0=',num2str(Parm{4}(1)),'\newline n_0/n_c=',num2str(Parm{4}(8)),...
    '\newline \omega_0=',num2str(w0),'\newline \lambda_0=',num2str(Data{2}(15)),...
    '\newline nw=',num2str(Parm{4}(16)),'\newline pulse shape=',num2str(Parm{4}(17)),...
    '\newline ndg=',num2str(ndg0),'\newline ndt=',num2str(ndt0),...
    '\newline np=',num2str(np0),'\newline lg=',num2str(Data{2}(3)), ' (',num2str(Data{2}(3)/Data{2}(15)),' \lambda _0)',...
    '\newline thickness=',num2str(Parm{4}(21)),' \lambda_0','\newline tmax=',num2str(Data{2}(4)),...
    '\newline xgleft=',num2str(Data{2}(9)),'  xgright=',num2str(Data{2}(10)),...
    '\newline \Delta g=',num2str(deltag),'  \Delta t=',num2str(deltat),...
    '\newline sample ndg=',num2str(sample_ndg),'  sample inc g=',num2str(sample_inc_g),...
    '\newline sample ndt=',num2str(sample_ndt),'  sample inc t=',num2str(sample_inc_t),...
    '\newline field seeing (transm. and ref.)=',num2str(Parm{4}(11)),' \lambda_0',...
    '\newline Fg MovingIons=',num2str(Parm{4}(7)),'  ionmass=',num2str(Parm{4}(12)),...
    '\newline sampled particles=',num2str(simulation_steps_p),...
    '\newline sampled spatial points=',num2str(simulation_steps_g),'\newline sampled temporal points=',num2str(simulation_steps_t)],...
    'HorizontalAlignment','center','BackgroundColor',[1 1 1],'FontSize',14);

figure
A=importdata('output3.txt', '\t');
%n=floor(ndg0/sample_div_ndg)%819%floor(size(A)/10);
n=simulation_steps_g;
%n=4
m=length(A(:,3))/n
length(A(:,3))
X=zeros(1,n);
Y=zeros(1,m);
for i=1:n
    X(i)=A(i,1);
end
for i=1:m
    Y(i)=A((i-1)*n+1,2);
end
X;
Y;
B=reshape(A(:,3),n,m);
length(X)
length(Y)
size(B)
iB = find (B<=-30);
B(iB) = NaN;
%surf(log10(abs(B)),'FaceColor','interp',...
surf(Y,X,B,'FaceColor','interp',...
'EdgeColor','none',...
 'FaceLighting','phong')
%daspect([5 5 1])
%axis([0 20 0 5 -12 0])
axis tight
colorbar
title('Density in  ne/nc')
xlabel('time (optical cycles)')
ylabel('spatial position')
view(2)
text(0.5,(max(Y))/2+4*Parm{4}(21),[' E_0=',num2str(Parm{4}(1)),'\newline n_0/n_c=',num2str(Parm{4}(8)),...
    '\newline ndg=',num2str(ndg0),'  sampled g=',num2str(simulation_steps_g),...
    '\newline ndt=',num2str(ndt0),'  sampled t=',num2str(simulation_steps_t),...
    '\newline np=',num2str(np0),'  sampled part=',num2str(simulation_steps_p),...
    '\newline \Delta g=',num2str(deltag),'  \Delta t=',num2str(deltat)],...
    'HorizontalAlignment','left','BackgroundColor',[1 .9 .9]);
%text(.1,.5,['\fontsize{16}black {\color{magenta}magenta '...
%'\color[rgb]{0 .5 .5}teal \color{red}red} black again'])
%text('Interpreter','latex',...
% 'String',['$$E_0=$$ ',num2str(Parm{4}(1)),'\\ $$n_0/n_c=$$ ',num2str(Parm{4}(8)),'\\ $$ndg=$$ ',num2str(ndg0)],...
% 'Position',[3 3],...
% 'FontSize',12,'BackgroundColor',[.9 .8 .8])

figure
Ac=importdata('output4.txt', '\t');
%nc=floor(ndg0/sample_div_ndg)%819%floor(size(A)/10);
nc=simulation_steps_g;
%n=4
mc=length(Ac(:,3))/nc
length(Ac(:,3))
Xc=zeros(1,nc);
Yc=zeros(1,mc);
for i=1:nc
    Xc(i)=Ac(i,1);
end
for i=1:mc
    Yc(i)=Ac((i-1)*nc+1,2);
end
Xc;
Yc;
Bc=reshape(Ac(:,3),nc,mc);
iBc = find (Bc<=-22);
Bc(iBc) = -22;
iBc = find (Bc>=22);
Bc(iBc) = 22;

length(Xc)
length(Yc)
size(Bc)
%surf(log10(abs(B)),'FaceColor','interp',...
surf(Yc,Xc,Bc,'FaceColor','interp',...
'EdgeColor','none',...
 'FaceLighting','phong')
%daspect([5 5 1])
axis tight
colorbar
title('Longitudinal field in u.a.')
xlabel('time (optical cycles)')
ylabel('spatial position')
view(2)
text(0.5,(max(Y))/2+4*Parm{4}(21),[' E_0=',num2str(Parm{4}(1)),'\newline n_0/n_c=',num2str(Parm{4}(8)),...
    '\newline ndg=',num2str(ndg0),'  sampled g=',num2str(simulation_steps_g),...
    '\newline ndt=',num2str(ndt0),'  sampled t=',num2str(simulation_steps_t),...
    '\newline np=',num2str(np0),'  sampled part=',num2str(simulation_steps_p),...
    '\newline \Delta g=',num2str(deltag),'  \Delta t=',num2str(deltat)],...
    'HorizontalAlignment','left','BackgroundColor',[.9 1 .9]);
%ylim([300 500])
%camlight lef




figure

fid2=fopen('output2.txt', 'r');

%Al leer el fichero me salto las primeras np0 filas, ya que es la
%distribucion incial y que engloba a todas las particulas.
vel=textscan(fid2, '%f\t%f\t%f\t%f', 'headerLines', np0);
vel{2}=vel{2}/lightC_au;
vel{3}=vel{3}/lightC_au;
vel{4}=vel{4}/lightC_au;
fclose(fid2);

maxvx=max(vel{2})
minvx=min(vel{2})
vsteps=300;
dvx=(maxvx-minvx)/vsteps
VX=zeros(simulation_steps_t+1, vsteps+1);
for i=1:simulation_steps_t*simulation_steps_p
    v_cell=(floor((vel{2}(i)+abs(minvx))/dvx))+1;
    VX(floor(i/simulation_steps_p)+1,v_cell)=VX(floor(i/simulation_steps_p)+1,v_cell)+1;
end
Xvx=zeros(1,simulation_steps_t+1);
for i=1:simulation_steps_t
    Xvx(i)=vel{1}((i-1)*simulation_steps_p+1);
end
Xvx(simulation_steps_t+1)=Xvx(simulation_steps_t)+sample_inc_t;
Yvx=minvx:dvx:maxvx;
length(Xvx)
length(Yvx)
size(VX)


iVX = find (VX>=100);
VX(iVX) = 100;


surf(Yvx,Xvx,VX,'FaceColor','interp',...
'EdgeColor','none',...
 'FaceLighting','phong')
axis tight
%zlim([0 100])
colorbar
title('Velocity distribution in the x-direction')
ylabel('time (optical cycles)')
xlabel('velocity (in c units)')
view(2)

%Distribucion inicial de velocidades
figure
fid2i=fopen('output2.txt', 'r');
Initvel=textscan(fid2i, '%f\t%f\t%f\t%f', np0);
Initvel{2}=Initvel{2}/lightC_au;
Initvel{3}=Initvel{3}/lightC_au;
Initvel{4}=Initvel{4}/lightC_au;
fclose(fid2i);
vIsteps=300;

subplot(2,2,1);
maxInitvx=max(Initvel{2})
minInitvx=min(Initvel{2})
vIsteps=300;
dIvx=(maxInitvx-minInitvx)/vIsteps
InitVx=zeros(1,vsteps+1);
for i=1:np0
    Initv_cell=(floor((Initvel{2}(i)+abs(minInitvx))/dIvx))+1;
    InitVx(Initv_cell)=InitVx(Initv_cell)+1;
end
YIvx=minInitvx:dIvx:maxInitvx;
length(YIvx);
length(InitVx);
plot(YIvx, InitVx, '-b')
title('Initial velocity distribution in the x-direction')
ylabel('counts')
xlabel('velocity (in c units)')

subplot(2,2,2);
maxInitvy=max(Initvel{3})
minInitvy=min(Initvel{3})
dIvy=(maxInitvy-minInitvy)/vIsteps
InitVy=zeros(1,vsteps+1);
for i=1:np0
    Initv_cell=(floor((Initvel{3}(i)+abs(minInitvy))/dIvy))+1;
    InitVy(Initv_cell)=InitVy(Initv_cell)+1;
end
YIvy=minInitvy:dIvy:maxInitvy;
length(YIvy);
length(InitVy);
plot(YIvy, InitVy, '-r')
title('Initial velocity distribution in the y-direction')
ylabel('counts')
xlabel('velocity (in c units)')

subplot(2,2,3);
maxInitvz=max(Initvel{4})
minInitvz=min(Initvel{4})
dIvz=(maxInitvz-minInitvz)/vIsteps
InitVz=zeros(1,vsteps+1);
for i=1:np0
    Initv_cell=(floor((Initvel{4}(i)+abs(minInitvz))/dIvz))+1;
    InitVz(Initv_cell)=InitVz(Initv_cell)+1;
end
YIvz=minInitvz:dIvz:maxInitvz;
length(YIvz);
length(InitVz);
plot(YIvz, InitVz, '-g')
title('Initial velocity distribution in the z-direction')
ylabel('counts')
xlabel('velocity (in c units)')




%listplot(so2{1}, so2{3},'LineWidth',1)
%xlabel('Time steps')
%ylabel('Position')
%title('Trayectories of the particles in time')


%figure
%listplot(so2{3}, so2{4},'LineWidth',1)
%xlabel('Posotion x')
%ylabel('Velocity x')
%title('Position vs velocity in the x-axis')

%fid3=fopen('datapic4.txt', 'r');

%pos=textscan(fid3, '%f\t%f');
%plot(pos{1}, pos{2},'Linestyle','none','Marker','.','MarkerSize',1);
%figure
%plot(pos{1}, pos{2});

%axis([0 4*1024 -0.5 0.5])
%xlabel('Time')
%ylabel('Position')
%title('Trayectories of the particles in time')


%fclose(fid3);